{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Anisotropic models\n",
    "\n",
    "In this tutorial, we demonstrate how to perform estimation with anisotropic variograms using GeoStats.jl.\n",
    "\n",
    "Before we proceed, please install the following packages:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\u001b[1m\u001b[36mINFO: \u001b[39m\u001b[22m\u001b[36mPackage GeoStats is already installed\n",
      "\u001b[39m\u001b[1m\u001b[36mINFO: \u001b[39m\u001b[22m\u001b[36mPackage Plots is already installed\n",
      "\u001b[39m\u001b[1m\u001b[36mINFO: \u001b[39m\u001b[22m\u001b[36mPackage GR is already installed\n",
      "\u001b[39m"
     ]
    }
   ],
   "source": [
    "for pkg in [\"GeoStats\", \"Plots\", \"GR\"]\n",
    "    Pkg.add(pkg)\n",
    "end\n",
    "\n",
    "# make sure this tutorial is reproducible\n",
    "srand(2000);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Ellipsoid distance\n",
    "\n",
    "Anisotropy can be thought of as a deformation of space with an ellipsoid distance. The semiaxes of the ellipsoid determine the preferential directions of the field and their lengths characterize the anisotropy ratio. In GeoStats.jl, all variogram models (empirical and theoretical) support a custom distance function that\n",
    "can be used to model anisotropy.\n",
    "\n",
    "A variogram object $\\gamma$ can be evaluated as an isotropic model $\\gamma(h)$ or as a (possibly) anisotropic model $\\gamma(x,y)$. For the Euclidean distance (the default), these two operations match $\\gamma(x,y) = \\gamma(h)$ in all directions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "true"
      ],
      "text/plain": [
       "true"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using GeoStats\n",
    "\n",
    "γ = GaussianVariogram()\n",
    "\n",
    "γ([1.,0.], [0.,0.]) ≈ γ(1.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If instead of an Euclidean ball, we use an ellipsoid with different semiaxes, the operation $\\gamma(x,y)$ becomes a function of the direction $x-y$. For example, we can create an ellipsoid distance aligned with the coordinate system where the major semiaxis has twice the size of the minor semiaxis:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "true"
      ],
      "text/plain": [
       "true"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using Distances\n",
    "\n",
    "d = Ellipsoidal([2.,1.],[0.])\n",
    "\n",
    "γaniso = GaussianVariogram(distance=d)\n",
    "\n",
    "γaniso([1.,0.],[0.,0.]) ≠ γaniso([0.,1.],[0.,0.])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Effects on estimation\n",
    "\n",
    "Now that we know how to construct anisotropic variograms, we can investigate the effect of varying the anisotropy ratio and alignement angle on estimation results. We start by generating some random data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<?xml version=\"1.0\" encoding=\"utf-8\"?>\n",
       "<svg xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\" width=\"600\" height=\"400\" viewBox=\"0 0 600 400\">\n",
       "<defs>\n",
       "  <clipPath id=\"clip00\">\n",
       "    <rect x=\"0\" y=\"0\" width=\"600\" height=\"400\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<polygon clip-path=\"url(#clip8000)\" points=\"\n",
       "0,400 600,400 600,0 0,0 \n",
       "  \" fill=\"#ffffff\" fill-opacity=\"1\"/>\n",
       "<defs>\n",
       "  <clipPath id=\"clip01\">\n",
       "    <rect x=\"120\" y=\"0\" width=\"421\" height=\"400\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<polygon clip-path=\"url(#clip8000)\" points=\"\n",
       "35.5169,375.813 580.315,375.813 580.315,11.811 35.5169,11.811 \n",
       "  \" fill=\"#ffffff\" fill-opacity=\"1\"/>\n",
       "<defs>\n",
       "  <clipPath id=\"clip02\">\n",
       "    <rect x=\"35\" y=\"11\" width=\"546\" height=\"365\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<polyline clip-path=\"url(#clip8002)\" style=\"stroke:#000000; stroke-width:0.5; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  138.921,370.353 138.921,17.2711 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip8002)\" style=\"stroke:#000000; stroke-width:0.5; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  250.184,370.353 250.184,17.2711 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip8002)\" style=\"stroke:#000000; stroke-width:0.5; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  361.447,370.353 361.447,17.2711 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip8002)\" style=\"stroke:#000000; stroke-width:0.5; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  472.71,370.353 472.71,17.2711 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip8002)\" style=\"stroke:#000000; stroke-width:0.5; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  43.6889,304.886 572.143,304.886 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip8002)\" style=\"stroke:#000000; stroke-width:0.5; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  43.6889,233.772 572.143,233.772 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip8002)\" style=\"stroke:#000000; stroke-width:0.5; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  43.6889,162.658 572.143,162.658 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip8002)\" style=\"stroke:#000000; stroke-width:0.5; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  43.6889,91.5434 572.143,91.5434 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip8002)\" style=\"stroke:#000000; stroke-width:0.5; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  43.6889,20.4291 572.143,20.4291 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip8000)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  35.5169,375.813 580.315,375.813 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip8000)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  35.5169,375.813 35.5169,11.811 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip8000)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  138.921,375.813 138.921,370.353 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip8000)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  250.184,375.813 250.184,370.353 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip8000)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  361.447,375.813 361.447,370.353 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip8000)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  472.71,375.813 472.71,370.353 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip8000)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  35.5169,304.886 43.6889,304.886 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip8000)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  35.5169,233.772 43.6889,233.772 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip8000)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  35.5169,162.658 43.6889,162.658 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip8000)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  35.5169,91.5434 43.6889,91.5434 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip8000)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  35.5169,20.4291 43.6889,20.4291 \n",
       "  \"/>\n",
       "<g clip-path=\"url(#clip8000)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:middle;\" transform=\"rotate(0, 138.921, 389.613)\" x=\"138.921\" y=\"389.613\">20</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip8000)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:middle;\" transform=\"rotate(0, 250.184, 389.613)\" x=\"250.184\" y=\"389.613\">40</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip8000)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:middle;\" transform=\"rotate(0, 361.447, 389.613)\" x=\"361.447\" y=\"389.613\">60</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip8000)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:middle;\" transform=\"rotate(0, 472.71, 389.613)\" x=\"472.71\" y=\"389.613\">80</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip8000)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:end;\" transform=\"rotate(0, 29.5169, 309.386)\" x=\"29.5169\" y=\"309.386\">20</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip8000)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:end;\" transform=\"rotate(0, 29.5169, 238.272)\" x=\"29.5169\" y=\"238.272\">40</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip8000)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:end;\" transform=\"rotate(0, 29.5169, 167.158)\" x=\"29.5169\" y=\"167.158\">60</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip8000)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:end;\" transform=\"rotate(0, 29.5169, 96.0434)\" x=\"29.5169\" y=\"96.0434\">80</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip8000)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:end;\" transform=\"rotate(0, 29.5169, 24.9291)\" x=\"29.5169\" y=\"24.9291\">100</text>\n",
       "</g>\n",
       "<circle clip-path=\"url(#clip8002)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"356.238\" cy=\"166.735\" r=\"4\"/>\n",
       "<circle clip-path=\"url(#clip8002)\" style=\"fill:#fafc9d; stroke:none; fill-opacity:1\" cx=\"356.238\" cy=\"166.735\" r=\"3\"/>\n",
       "<circle clip-path=\"url(#clip8002)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"103.66\" cy=\"200.519\" r=\"4\"/>\n",
       "<circle clip-path=\"url(#clip8002)\" style=\"fill:#f7860e; stroke:none; fill-opacity:1\" cx=\"103.66\" cy=\"200.519\" r=\"3\"/>\n",
       "<circle clip-path=\"url(#clip8002)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"327.721\" cy=\"194.135\" r=\"4\"/>\n",
       "<circle clip-path=\"url(#clip8002)\" style=\"fill:#f99008; stroke:none; fill-opacity:1\" cx=\"327.721\" cy=\"194.135\" r=\"3\"/>\n",
       "<circle clip-path=\"url(#clip8002)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"448.365\" cy=\"157.438\" r=\"4\"/>\n",
       "<circle clip-path=\"url(#clip8002)\" style=\"fill:#be3852; stroke:none; fill-opacity:1\" cx=\"448.365\" cy=\"157.438\" r=\"3\"/>\n",
       "<circle clip-path=\"url(#clip8002)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"439.207\" cy=\"204.281\" r=\"4\"/>\n",
       "<circle clip-path=\"url(#clip8002)\" style=\"fill:#b63457; stroke:none; fill-opacity:1\" cx=\"439.207\" cy=\"204.281\" r=\"3\"/>\n",
       "<circle clip-path=\"url(#clip8002)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"56.556\" cy=\"279.059\" r=\"4\"/>\n",
       "<circle clip-path=\"url(#clip8002)\" style=\"fill:#8e2368; stroke:none; fill-opacity:1\" cx=\"56.556\" cy=\"279.059\" r=\"3\"/>\n",
       "<circle clip-path=\"url(#clip8002)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"541.215\" cy=\"108.799\" r=\"4\"/>\n",
       "<circle clip-path=\"url(#clip8002)\" style=\"fill:#d84c3d; stroke:none; fill-opacity:1\" cx=\"541.215\" cy=\"108.799\" r=\"3\"/>\n",
       "<circle clip-path=\"url(#clip8002)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"419.231\" cy=\"62.0243\" r=\"4\"/>\n",
       "<circle clip-path=\"url(#clip8002)\" style=\"fill:#c73e4b; stroke:none; fill-opacity:1\" cx=\"419.231\" cy=\"62.0243\" r=\"3\"/>\n",
       "<circle clip-path=\"url(#clip8002)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"56.4461\" cy=\"233.185\" r=\"4\"/>\n",
       "<circle clip-path=\"url(#clip8002)\" style=\"fill:#000003; stroke:none; fill-opacity:1\" cx=\"56.4461\" cy=\"233.185\" r=\"3\"/>\n",
       "<circle clip-path=\"url(#clip8002)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"50.9357\" cy=\"85.1541\" r=\"4\"/>\n",
       "<circle clip-path=\"url(#clip8002)\" style=\"fill:#1d0b43; stroke:none; fill-opacity:1\" cx=\"50.9357\" cy=\"85.1541\" r=\"3\"/>\n",
       "<circle clip-path=\"url(#clip8002)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"376.599\" cy=\"177.583\" r=\"4\"/>\n",
       "<circle clip-path=\"url(#clip8002)\" style=\"fill:#450a68; stroke:none; fill-opacity:1\" cx=\"376.599\" cy=\"177.583\" r=\"3\"/>\n",
       "<circle clip-path=\"url(#clip8002)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"308.019\" cy=\"62.2716\" r=\"4\"/>\n",
       "<circle clip-path=\"url(#clip8002)\" style=\"fill:#972765; stroke:none; fill-opacity:1\" cx=\"308.019\" cy=\"62.2716\" r=\"3\"/>\n",
       "<circle clip-path=\"url(#clip8002)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"197.048\" cy=\"52.8341\" r=\"4\"/>\n",
       "<circle clip-path=\"url(#clip8002)\" style=\"fill:#eb6527; stroke:none; fill-opacity:1\" cx=\"197.048\" cy=\"52.8341\" r=\"3\"/>\n",
       "<circle clip-path=\"url(#clip8002)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"155.196\" cy=\"61.6236\" r=\"4\"/>\n",
       "<circle clip-path=\"url(#clip8002)\" style=\"fill:#f5d645; stroke:none; fill-opacity:1\" cx=\"155.196\" cy=\"61.6236\" r=\"3\"/>\n",
       "<circle clip-path=\"url(#clip8002)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"564.896\" cy=\"232.04\" r=\"4\"/>\n",
       "<circle clip-path=\"url(#clip8002)\" style=\"fill:#721a6d; stroke:none; fill-opacity:1\" cx=\"564.896\" cy=\"232.04\" r=\"3\"/>\n",
       "<circle clip-path=\"url(#clip8002)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"448.497\" cy=\"35.7554\" r=\"4\"/>\n",
       "<circle clip-path=\"url(#clip8002)\" style=\"fill:#f9fc9d; stroke:none; fill-opacity:1\" cx=\"448.497\" cy=\"35.7554\" r=\"3\"/>\n",
       "<circle clip-path=\"url(#clip8002)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"535.301\" cy=\"237.047\" r=\"4\"/>\n",
       "<circle clip-path=\"url(#clip8002)\" style=\"fill:#f6d23f; stroke:none; fill-opacity:1\" cx=\"535.301\" cy=\"237.047\" r=\"3\"/>\n",
       "<circle clip-path=\"url(#clip8002)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"60.531\" cy=\"215.202\" r=\"4\"/>\n",
       "<circle clip-path=\"url(#clip8002)\" style=\"fill:#c53d4d; stroke:none; fill-opacity:1\" cx=\"60.531\" cy=\"215.202\" r=\"3\"/>\n",
       "<circle clip-path=\"url(#clip8002)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"470.896\" cy=\"312.015\" r=\"4\"/>\n",
       "<circle clip-path=\"url(#clip8002)\" style=\"fill:#de5337; stroke:none; fill-opacity:1\" cx=\"470.896\" cy=\"312.015\" r=\"3\"/>\n",
       "<circle clip-path=\"url(#clip8002)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"432.397\" cy=\"336.997\" r=\"4\"/>\n",
       "<circle clip-path=\"url(#clip8002)\" style=\"fill:#fbb317; stroke:none; fill-opacity:1\" cx=\"432.397\" cy=\"336.997\" r=\"3\"/>\n",
       "<circle clip-path=\"url(#clip8002)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"436.869\" cy=\"39.7225\" r=\"4\"/>\n",
       "<circle clip-path=\"url(#clip8002)\" style=\"fill:#0d0828; stroke:none; fill-opacity:1\" cx=\"436.869\" cy=\"39.7225\" r=\"3\"/>\n",
       "<circle clip-path=\"url(#clip8002)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"318.498\" cy=\"24.3023\" r=\"4\"/>\n",
       "<circle clip-path=\"url(#clip8002)\" style=\"fill:#1a0b3f; stroke:none; fill-opacity:1\" cx=\"318.498\" cy=\"24.3023\" r=\"3\"/>\n",
       "<circle clip-path=\"url(#clip8002)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"116.627\" cy=\"136.864\" r=\"4\"/>\n",
       "<circle clip-path=\"url(#clip8002)\" style=\"fill:#a42c60; stroke:none; fill-opacity:1\" cx=\"116.627\" cy=\"136.864\" r=\"3\"/>\n",
       "<circle clip-path=\"url(#clip8002)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"359.66\" cy=\"22.113\" r=\"4\"/>\n",
       "<circle clip-path=\"url(#clip8002)\" style=\"fill:#fba40b; stroke:none; fill-opacity:1\" cx=\"359.66\" cy=\"22.113\" r=\"3\"/>\n",
       "<circle clip-path=\"url(#clip8002)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"436.69\" cy=\"234.639\" r=\"4\"/>\n",
       "<circle clip-path=\"url(#clip8002)\" style=\"fill:#f2f380; stroke:none; fill-opacity:1\" cx=\"436.69\" cy=\"234.639\" r=\"3\"/>\n",
       "<circle clip-path=\"url(#clip8002)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"317.358\" cy=\"165.123\" r=\"4\"/>\n",
       "<circle clip-path=\"url(#clip8002)\" style=\"fill:#5c116d; stroke:none; fill-opacity:1\" cx=\"317.358\" cy=\"165.123\" r=\"3\"/>\n",
       "<circle clip-path=\"url(#clip8002)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"109.729\" cy=\"56.373\" r=\"4\"/>\n",
       "<circle clip-path=\"url(#clip8002)\" style=\"fill:#e55c2f; stroke:none; fill-opacity:1\" cx=\"109.729\" cy=\"56.373\" r=\"3\"/>\n",
       "<circle clip-path=\"url(#clip8002)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"122.563\" cy=\"177.724\" r=\"4\"/>\n",
       "<circle clip-path=\"url(#clip8002)\" style=\"fill:#2b0a55; stroke:none; fill-opacity:1\" cx=\"122.563\" cy=\"177.724\" r=\"3\"/>\n",
       "<circle clip-path=\"url(#clip8002)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"445.286\" cy=\"54.4649\" r=\"4\"/>\n",
       "<circle clip-path=\"url(#clip8002)\" style=\"fill:#fac027; stroke:none; fill-opacity:1\" cx=\"445.286\" cy=\"54.4649\" r=\"3\"/>\n",
       "<circle clip-path=\"url(#clip8002)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"166.821\" cy=\"280.052\" r=\"4\"/>\n",
       "<circle clip-path=\"url(#clip8002)\" style=\"fill:#f7840f; stroke:none; fill-opacity:1\" cx=\"166.821\" cy=\"280.052\" r=\"3\"/>\n",
       "<circle clip-path=\"url(#clip8002)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"362.804\" cy=\"365.511\" r=\"4\"/>\n",
       "<circle clip-path=\"url(#clip8002)\" style=\"fill:#fcfea4; stroke:none; fill-opacity:1\" cx=\"362.804\" cy=\"365.511\" r=\"3\"/>\n",
       "<circle clip-path=\"url(#clip8002)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"79.3315\" cy=\"170.375\" r=\"4\"/>\n",
       "<circle clip-path=\"url(#clip8002)\" style=\"fill:#1e0b45; stroke:none; fill-opacity:1\" cx=\"79.3315\" cy=\"170.375\" r=\"3\"/>\n",
       "<circle clip-path=\"url(#clip8002)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"472.57\" cy=\"234.362\" r=\"4\"/>\n",
       "<circle clip-path=\"url(#clip8002)\" style=\"fill:#831f6b; stroke:none; fill-opacity:1\" cx=\"472.57\" cy=\"234.362\" r=\"3\"/>\n",
       "<circle clip-path=\"url(#clip8002)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"190.436\" cy=\"156.982\" r=\"4\"/>\n",
       "<circle clip-path=\"url(#clip8002)\" style=\"fill:#260b51; stroke:none; fill-opacity:1\" cx=\"190.436\" cy=\"156.982\" r=\"3\"/>\n",
       "<circle clip-path=\"url(#clip8002)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"58.8945\" cy=\"145.808\" r=\"4\"/>\n",
       "<circle clip-path=\"url(#clip8002)\" style=\"fill:#f7f994; stroke:none; fill-opacity:1\" cx=\"58.8945\" cy=\"145.808\" r=\"3\"/>\n",
       "<circle clip-path=\"url(#clip8002)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"328.228\" cy=\"262.107\" r=\"4\"/>\n",
       "<circle clip-path=\"url(#clip8002)\" style=\"fill:#9a2864; stroke:none; fill-opacity:1\" cx=\"328.228\" cy=\"262.107\" r=\"3\"/>\n",
       "<circle clip-path=\"url(#clip8002)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"316.994\" cy=\"63.4827\" r=\"4\"/>\n",
       "<circle clip-path=\"url(#clip8002)\" style=\"fill:#1a0b40; stroke:none; fill-opacity:1\" cx=\"316.994\" cy=\"63.4827\" r=\"3\"/>\n",
       "<circle clip-path=\"url(#clip8002)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"117.794\" cy=\"129.653\" r=\"4\"/>\n",
       "<circle clip-path=\"url(#clip8002)\" style=\"fill:#02020e; stroke:none; fill-opacity:1\" cx=\"117.794\" cy=\"129.653\" r=\"3\"/>\n",
       "<circle clip-path=\"url(#clip8002)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"79.8438\" cy=\"298.747\" r=\"4\"/>\n",
       "<circle clip-path=\"url(#clip8002)\" style=\"fill:#c63d4c; stroke:none; fill-opacity:1\" cx=\"79.8438\" cy=\"298.747\" r=\"3\"/>\n",
       "<circle clip-path=\"url(#clip8002)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"452.367\" cy=\"87.4041\" r=\"4\"/>\n",
       "<circle clip-path=\"url(#clip8002)\" style=\"fill:#e55c2f; stroke:none; fill-opacity:1\" cx=\"452.367\" cy=\"87.4041\" r=\"3\"/>\n",
       "<circle clip-path=\"url(#clip8002)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"290.795\" cy=\"93.2022\" r=\"4\"/>\n",
       "<circle clip-path=\"url(#clip8002)\" style=\"fill:#fba40b; stroke:none; fill-opacity:1\" cx=\"290.795\" cy=\"93.2022\" r=\"3\"/>\n",
       "<circle clip-path=\"url(#clip8002)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"105.37\" cy=\"161.231\" r=\"4\"/>\n",
       "<circle clip-path=\"url(#clip8002)\" style=\"fill:#ef6e20; stroke:none; fill-opacity:1\" cx=\"105.37\" cy=\"161.231\" r=\"3\"/>\n",
       "<circle clip-path=\"url(#clip8002)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"293.502\" cy=\"164.109\" r=\"4\"/>\n",
       "<circle clip-path=\"url(#clip8002)\" style=\"fill:#c53d4d; stroke:none; fill-opacity:1\" cx=\"293.502\" cy=\"164.109\" r=\"3\"/>\n",
       "<circle clip-path=\"url(#clip8002)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"423.647\" cy=\"340.492\" r=\"4\"/>\n",
       "<circle clip-path=\"url(#clip8002)\" style=\"fill:#520e6c; stroke:none; fill-opacity:1\" cx=\"423.647\" cy=\"340.492\" r=\"3\"/>\n",
       "<circle clip-path=\"url(#clip8002)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"97.9806\" cy=\"200.316\" r=\"4\"/>\n",
       "<circle clip-path=\"url(#clip8002)\" style=\"fill:#7f1e6c; stroke:none; fill-opacity:1\" cx=\"97.9806\" cy=\"200.316\" r=\"3\"/>\n",
       "<circle clip-path=\"url(#clip8002)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"94.9799\" cy=\"361.234\" r=\"4\"/>\n",
       "<circle clip-path=\"url(#clip8002)\" style=\"fill:#e96329; stroke:none; fill-opacity:1\" cx=\"94.9799\" cy=\"361.234\" r=\"3\"/>\n",
       "<circle clip-path=\"url(#clip8002)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"72.6234\" cy=\"65.5422\" r=\"4\"/>\n",
       "<circle clip-path=\"url(#clip8002)\" style=\"fill:#ea6428; stroke:none; fill-opacity:1\" cx=\"72.6234\" cy=\"65.5422\" r=\"3\"/>\n",
       "<circle clip-path=\"url(#clip8002)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"516.172\" cy=\"228.082\" r=\"4\"/>\n",
       "<circle clip-path=\"url(#clip8002)\" style=\"fill:#f78310; stroke:none; fill-opacity:1\" cx=\"516.172\" cy=\"228.082\" r=\"3\"/>\n",
       "<circle clip-path=\"url(#clip8002)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"67.7617\" cy=\"193.353\" r=\"4\"/>\n",
       "<circle clip-path=\"url(#clip8002)\" style=\"fill:#b0315a; stroke:none; fill-opacity:1\" cx=\"67.7617\" cy=\"193.353\" r=\"3\"/>\n",
       "<circle clip-path=\"url(#clip8002)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"238.321\" cy=\"231.205\" r=\"4\"/>\n",
       "<circle clip-path=\"url(#clip8002)\" style=\"fill:#360961; stroke:none; fill-opacity:1\" cx=\"238.321\" cy=\"231.205\" r=\"3\"/>\n",
       "<polygon clip-path=\"url(#clip8000)\" points=\"\n",
       "479.251,62.931 562.315,62.931 562.315,32.691 479.251,32.691 \n",
       "  \" fill=\"#ffffff\" fill-opacity=\"1\"/>\n",
       "<polyline clip-path=\"url(#clip8000)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  479.251,62.931 562.315,62.931 562.315,32.691 479.251,32.691 479.251,62.931 \n",
       "  \"/>\n",
       "<circle clip-path=\"url(#clip8000)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"506.251\" cy=\"47.811\" r=\"6\"/>\n",
       "<circle clip-path=\"url(#clip8000)\" style=\"fill:#fafc9d; stroke:none; fill-opacity:1\" cx=\"506.251\" cy=\"47.811\" r=\"5\"/>\n",
       "<g clip-path=\"url(#clip8000)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:start;\" transform=\"rotate(0, 527.251, 52.311)\" x=\"527.251\" y=\"52.311\">data</text>\n",
       "</g>\n",
       "</svg>\n"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using Plots; gr()\n",
    "\n",
    "dim, nobs = 2, 50\n",
    "\n",
    "X = 100*rand(dim, nobs)\n",
    "z = rand(nobs)\n",
    "\n",
    "scatter(X[1,:], X[2,:], c=z, label=\"data\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and by defining an estimation problem:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2D EstimationProblem\n",
       "  data:      50×3 GeoDataFrame (x and y)\n",
       "  domain:    100×100 RegularGrid{Float64,2}\n",
       "  variables: variable (Float64)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "geodata = GeoDataFrame(DataFrames.DataFrame(x=X[1,:],y=X[2,:],variable=z), [:x,:y])\n",
    "domain = RegularGrid{Float64}(100,100)\n",
    "problem = EstimationProblem(geodata, domain, :variable)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, we vary the anisotropy ratio with an ellipsoid that is aligned with the coordinate system:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\u001b[1m\u001b[36mINFO: \u001b[39m\u001b[22m\u001b[36mSaved animation to /home/juliohm/.julia/v0.6/GeoStats/examples/figs/anisotropy_ratio.gif\n",
      "\u001b[39m"
     ]
    },
    {
     "data": {
      "text/html": [
       "<img src=\"figs/anisotropy_ratio.gif?0.943356905739213>\" />"
      ],
      "text/plain": [
       "Plots.AnimatedGif(\"/home/juliohm/.julia/v0.6/GeoStats/examples/figs/anisotropy_ratio.gif\")"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "anim = @animate for a=linspace(1,10,10)\n",
    "    d = Ellipsoidal([a,1.], [0.])\n",
    "    \n",
    "    γ = GaussianVariogram(range=5., distance=d)\n",
    "    \n",
    "    solution = solve(problem, Kriging(:variable => @NT(variogram=γ)))\n",
    "    \n",
    "    plot(solution, size=(1000,400))\n",
    "end\n",
    "gif(anim, \"figs/anisotropy_ratio.gif\", fps=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Second, we fix the anisotropy ratio and vary the alignment angle:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\u001b[1m\u001b[36mINFO: \u001b[39m\u001b[22m\u001b[36mSaved animation to /home/juliohm/.julia/v0.6/GeoStats/examples/figs/anisotropy_angle.gif\n",
      "\u001b[39m"
     ]
    },
    {
     "data": {
      "text/html": [
       "<img src=\"figs/anisotropy_angle.gif?0.57557066645746>\" />"
      ],
      "text/plain": [
       "Plots.AnimatedGif(\"/home/juliohm/.julia/v0.6/GeoStats/examples/figs/anisotropy_angle.gif\")"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "anim = @animate for θ=linspace(0,2π,10)\n",
    "    d = Ellipsoidal([10.,1.], [θ])\n",
    "    \n",
    "    γ = GaussianVariogram(range=5., distance=d)\n",
    "    \n",
    "    solution = solve(problem, Kriging(:variable => @NT(variogram=γ)))\n",
    "    \n",
    "    plot(solution, size=(1000,400))\n",
    "end\n",
    "gif(anim, \"figs/anisotropy_angle.gif\", fps=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This experiment can be extended to 3D with the only difference being that ellipsoids therein are defined by 3 semiaxes and 3 angles. For example, the Euclidean distance in 3D can be recovered with a degenerated ellipsoid with equal semiaxes (i.e. sphere):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "true"
      ],
      "text/plain": [
       "true"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ellipsoid = Ellipsoidal([1.,1.,1.],[0.,0.,0.])\n",
    "euclidean = Euclidean()\n",
    "\n",
    "a, b = rand(3), rand(3)\n",
    "\n",
    "evaluate(ellipsoid, a, b) ≈ evaluate(euclidean, a, b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Please consult [Distances.jl](https://github.com/JuliaStats/Distances.jl) for more distance functions."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 0.6.1-pre",
   "language": "julia",
   "name": "julia-0.6"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
